Enunciado

En el archivo sat.txt se tienen datos de valores de espectros de pixels en una imagen de satélite utilizados para la predicción de la clase de suelo. Realice un análisis de los datos utilizando diferentes métodos y compárelos mediante Validación Cruzada. Luego aplíquelos a la muestra de test sat.tst y analice los resultados obtenidos.

Descripción de los Datos


La base de datos contiene los valores multi-espectrales de píxeles en regiones de 3x3 en una imágen satelital, y la clasificación asociada con el píxel central de cada región. El objetivo es predecir esta clasificación dados los valores multi-espectrales.

Esta base de datos se generó a partir del escáner multi-espectral LANDSAT. Un marco de imágenes Landsat MSS consta de cuatro imágenes digitales de la misma escena en diferentes bandas espectrales. Dos de estos están en la región visible (correspondiente aproximadamente a las regiones verdes y rojas del espectro visible) y dos están en el infrarrojo. Cada píxel es una palabra binaria de 8 bits, con 0 correspondiente al negro y 255 a blanco. La resolución espacial de un píxel es de unos 80 m x 80m. Cada imagen contiene 2340 x 3380 de tales píxeles.

La base de datos es una subárea (pequeña) de una escena, que consta de 82 x 100 píxeles. Cada línea de datos corresponde a una región de píxeles cuadrados 3x3 completamente contenida dentro de la subárea 82x100. Cada línea contiene los valores de píxeles en las cuatro bandas espectrales (convertidas en ASCII) de cada uno de los 9 píxeles en la región 3x3 y un número que indica la etiqueta de clasificación del píxel central.

La clase de cada píxel se codifica como un número, siendo estos:

En cada línea de datos, los cuatro valores espectrales para el píxel superior izquierdo se dan primero, seguidos por los cuatro valores espectrales para el píxel medio superior y luego los del píxel superior derecho, y así sucesivamente de modo que los píxeles se leen en secuencia de izquierda a derecha y de arriba a abajo. Por lo tanto, los cuatro valores espectrales para los píxeles centrales están dados por los atributos 17,18,19 y 20. Si lo desea, puede usar solo estos cuatro atributos, ignorando los demás.

Resolución

Obtención del Set de Datos

Dado que el enunciado menciona que pueden utilizarse solo los datos espectrales correspondientes al píxel central, buscamos extraerlos, junto con su clasificación, a una sola dataframe.

file_trn = "SAT_trn.txt"
data_trn <- read.table(file_trn, sep = "", header=F)
sat_trn_df <- data_trn      %>% 
  dplyr::select(17:20,37)  %>% 
  rename(B1=V17,
         B2=V18,
         B3=V19,
         B4=V20,
         C=V37)
file_tst = "SAT_tst.txt"
data_tst <- read.table(file_tst, sep = "", header=F)
sat_tst_df <- data_tst      %>% 
  dplyr::select(17:20,37)  %>% 
  rename(B1=V17,
         B2=V18,
         B3=V19,
         B4=V20,
         C=V37)
sat_trn_df$C[sat_trn_df$C == 7] <- 6
sat_tst_df$C[sat_tst_df$C == 7] <- 6
sat_trn_df_numeric <- data.frame(sat_trn_df)
sat_tst_df_numeric <- data.frame(sat_tst_df)
class_names <- c("rojo", "algodón", "gris","gris húmedo", "vegetación", "gris muy húmedo")
class_values <- c(1,2,3,4,5,6)
class_map <- setNames(class_names, class_values)
sat_trn_df$C <- class_map[sat_trn_df$C]
sat_tst_df$C <- class_map[sat_tst_df$C]

Visualización de los Datos

Contenido del Set de Datos

B1 B2 B3 B4 C
92 112 118 85 gris
84 103 104 81 gris
84 99 104 78 gris
84 99 104 81 gris
76 99 104 81 gris
76 99 108 85 gris
80 112 118 88 gris
80 107 113 85 gris
76 91 104 74 gris húmedo
76 95 100 78 gris húmedo

Podemos comprobar que tenemos valores acotados entre 0 y 255 para cuatro bandas espectrales y una categoría de suelo asignada en cada fila. Los contenidos pueden visualizarse mejor utilizando un gráfico de coordenadas paralelas.

Donde cada linea representa una observación para las 4 bandas. La normalización de los valores no es necesaria ya que \(B_1:B_4\) estan en la misma escala. Podemos ver que la distribución de suelo en los datos de entrenamiento y prueba es practicamente la misma.

Construcción de Modelos

sat_trn_df <- data.frame(sat_trn_df_numeric)
sat_tst_df <- data.frame(sat_tst_df_numeric)
sat_trn_df$C <- as.factor(sat_trn_df$C)
sat_tst_df$C <- as.factor(sat_tst_df$C)
set.seed(20)

Árbol de Decisión

Librería Tree

# Entrenamiento del Modelo
n <- nrow(sat_trn_df)
sat_tree <- tree(C ~ B1 + B2 + B3 + B4, 
                 data= sat_trn_df, 
                 control=tree.control(
                   nobs= nrow(sat_trn_df), 
                   mincut = 1, 
                   minsize = 2, 
                   mindev = 0)
                 )
# Resultados
sat_tree_summary <- summary(sat_tree)
EEnt <- sat_tree_summary$misclass[1] / n
sat_tree_summary
## 
## Classification tree:
## tree(formula = C ~ B1 + B2 + B3 + B4, data = sat_trn_df, control = tree.control(nobs = nrow(sat_trn_df), 
##     mincut = 1, minsize = 2, mindev = 0))
## Number of terminal nodes:  761 
## Residual mean deviance:  0.1723 = 633 / 3674 
## Misclassification error rate: 0.04014 = 178 / 4435
sat_tree_CV = cv.tree(sat_tree, FUN = prune.misclass)
size <- sat_tree_CV$size
dev <- sat_tree_CV$dev
sp <-ggplot(data=data.frame(size,dev), aes(x=size, y=dev)) +
  geom_line(linetype="solid", color="blue", size=1.2)+
  geom_point(color="red", size=3)+
  labs(x= "Tamaño del Árbol", y = "Error de CV",title="Poda de Complejidad de Costos")
sp + scale_x_continuous(trans='log10')

# Agrego código para obtener con R el mejor Tamaño
best_size <- which.min(sat_tree_CV$dev) 
best_size
## [1] 39
# PQ Harcodea el 13?
sat_tree_prune = prune.misclass(sat_tree, best = 13)
plot(sat_tree_prune)
text(sat_tree_prune, pretty = 0)

CVErr_tree <- sat_tree_CV$dev[39] / n
CVErr_tree
## [1] 0.1815107

Librería RPart

set.seed(20)
# Entrenamiento del Modelo con 10-fold CV  
sat_rpart <-  rpart(C ~ B1 + B2 + B3 + B4, data= sat_trn_df, method = "class", cp=-1,xval=10)
sat_rpart_summary <- printcp(sat_rpart)
## 
## Classification tree:
## rpart(formula = C ~ B1 + B2 + B3 + B4, data = sat_trn_df, method = "class", 
##     cp = -1, xval = 10)
## 
## Variables actually used in tree construction:
## [1] B1 B2 B3 B4
## 
## Root node error: 3363/4435 = 0.75829
## 
## n= 4435 
## 
##             CP nsplit rel error  xerror      xstd
## 1   2.6167e-01      0   1.00000 1.00000 0.0084779
## 2   2.5335e-01      1   0.73833 0.76896 0.0097636
## 3   1.2905e-01      2   0.48498 0.48498 0.0095487
## 4   6.3634e-02      3   0.35593 0.35742 0.0088020
## 5   1.7544e-02      4   0.29230 0.29379 0.0082400
## 6   1.4570e-02      5   0.27475 0.28189 0.0081181
## 7   6.6905e-03      7   0.24561 0.24234 0.0076694
## 8   5.9471e-03     10   0.22034 0.23253 0.0075467
## 9   5.6497e-03     11   0.21439 0.22629 0.0074660
## 10  4.7577e-03     13   0.20309 0.21885 0.0073673
## 11  4.6090e-03     14   0.19833 0.20934 0.0072364
## 12  3.2709e-03     16   0.18912 0.20636 0.0071945
## 13  1.6354e-03     17   0.18585 0.19774 0.0070698
## 14  1.3381e-03     19   0.18258 0.19625 0.0070479
## 15  1.1894e-03     21   0.17990 0.19536 0.0070346
## 16  1.0903e-03     22   0.17871 0.19477 0.0070258
## 17  8.9206e-04     25   0.17544 0.19715 0.0070610
## 18  6.9382e-04     31   0.17009 0.19715 0.0070610
## 19  5.9471e-04     34   0.16800 0.19715 0.0070610
## 20  4.9559e-04     48   0.15968 0.19685 0.0070567
## 21  4.4603e-04     51   0.15819 0.19387 0.0070124
## 22  3.9647e-04     55   0.15641 0.19328 0.0070035
## 23  2.9735e-04     58   0.15522 0.19328 0.0070035
## 24  1.4868e-04     61   0.15433 0.20012 0.0071047
## 25  8.4958e-05     63   0.15403 0.20101 0.0071176
## 26  7.4338e-05     70   0.15343 0.20071 0.0071133
## 27  0.0000e+00     74   0.15314 0.19952 0.0070960
## 28 -1.0000e+00    197   0.15314 0.19952 0.0070960
plotcp(sat_rpart) 

which.min(sat_rpart$cptable[,4])
## 22 
## 22
cp<-sat_rpart$cptable[which.min(sat_rpart$cptable[,"xerror"]),"CP"]
cp
## [1] 0.0003964714
CVErr_rpart = sat_rpart$frame[1,'dev'] / n * sat_rpart$cptable[which.min(sat_rpart$cptable[,"xerror"]),"xerror"]
CVErr_rpart
## [1] 0.1465614
poda_rpart<-prune(sat_rpart,cp=cp)
rpart.plot(poda_rpart, cex = 0.5, extra = 0)

Random Forest

set.seed(20)
sat_rf<- randomForest(C~., data= sat_trn_df, 
                     mtry=sqrt(ncol(sat_trn_df)-1), importance=TRUE)

importance(sat_rf)
##            1        2         3        4         5        6
## B1 227.87903 29.96135 362.16076 58.21353  50.21856 48.37936
## B2 140.86987 43.00905  53.51638 37.62347 112.03571 41.93234
## B3  30.04423 10.40657  24.77468 15.11574  28.56435 32.31841
## B4  57.76868 24.80617  34.26451 35.82601  52.16778 74.12139
##    MeanDecreaseAccuracy MeanDecreaseGini
## B1            186.14425        1186.7218
## B2            119.11363         914.1809
## B3             48.88156         425.2869
## B4             91.94379         853.1770
varImpPlot(sat_rf)

CVErr_rf <- unname(sat_rf$err.rate[nrow(sat_rf$err.rate),1])
CVErr_rf
## [1] 0.1492672

Vecinos Más Cercanos

set.seed(20)
sat_knn <- train(C ~ ., data = sat_trn_df, 
                   method = "knn",
                   preProcess = c("center","scale"),
                   trControl =  trainControl(method = "cv", number = 10),
                   tuneGrid = expand.grid(k = 1:100),
                   metric = "Accuracy")
results <- sat_knn$results
ggplot(results, aes(x = k, y = Accuracy)) + 
  geom_line(color = "steelblue", size = 1.2) + 
  xlab("Valor de K") + 
  ylab("Exactitud") + 
  ggtitle("Rendimiento para diferentes valores de K") +
  theme_bw() + 
  theme(plot.title = element_text(face = "bold", size = 14, hjust = 0.5),
        axis.title.x = element_text(face = "bold", size = 12),
        axis.title.y = element_text(face = "bold", size = 12),
        axis.text = element_text(size = 10),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank())

best_k <- results$k[which.max(results$Accuracy)]
best_k
## [1] 43
CVErr_KNN <- 1-unname(sat_knn$results[best_k,"Accuracy"])
CVErr_KNN
## [1] 0.1400425

Regresión Logística Multinomial

summary(sat_mlr)
## Call:
## nnet::multinom(formula = .outcome ~ ., data = dat, decay = param$decay)
## 
## Coefficients:
##   (Intercept)        B1          B2          B3         B4
## 2   -9.921154 0.5867093 -0.63317495  0.29943995 -0.1122932
## 3  -36.648978 0.7617752 -0.02673514  0.02795903 -0.2534971
## 4  -10.253548 0.5886442 -0.07999703 -0.02004337 -0.2854203
## 5    7.594031 0.4558313 -0.52173889  0.21832973 -0.2099831
## 6    4.599893 0.5678520 -0.18940106 -0.05394639 -0.2887664
## 
## Std. Errors:
##   (Intercept)         B1         B2         B3         B4
## 2   1.6977543 0.05040409 0.04477873 0.05284647 0.04508065
## 3   1.2832527 0.03391253 0.03743092 0.04631104 0.04434272
## 4   0.8182467 0.03437754 0.03467115 0.04409015 0.04040293
## 5   0.9108139 0.03161938 0.03489501 0.04077884 0.03601165
## 6   0.7841282 0.03491517 0.03359844 0.04276092 0.03778148
## 
## Residual Deviance: 3755.425 
## AIC: 3805.425
#sat_mlr <- nnet::multinom(C ~ ., data = sat_trn_df)
#summary(sat_rl)
# Predict the class labels for the training data
confusion_matrix <- confusionMatrix(predict(sat_mlr, newdata = sat_trn_df), sat_trn_df$C)

CVErr_mlr2 <- 1 - unname(confusion_matrix$overall['Accuracy'])
CVErr_mlr2
## [1] 0.1569335
predicted_mlr <- predict(sat_mlr, newdata = sat_trn_df)

CVErr_mlr <- mean(predicted_mlr != sat_trn_df$C)
CVErr_mlr
## [1] 0.1569335

Análisis Discriminante Lineal

set.seed(20)
sat_LDA <- stepclass(C ~ B1 + B2 + B3 + B4, data= sat_trn_df, method="lda", criterion = "CR", direction = "both", improvement <- 0.01, fold = 10)
##  `stepwise classification', using 10-fold cross-validated correctness rate of method lda'.
## 4435 observations of 4 variables in 6 classes; direction: both
## stop criterion: improvement less than 1%.
## correctness rate: 0.56054;  in: "B2";  variables (1): B2 
## correctness rate: 0.79776;  in: "B1";  variables (2): B2, B1 
## correctness rate: 0.82616;  in: "B3";  variables (3): B2, B1, B3 
## 
##  hr.elapsed min.elapsed sec.elapsed 
##        0.00        0.00        3.13
sat_LDA
## method      : lda 
## final model : C ~ B1 + B2 + B3
## <environment: 0x000001f8ba1e8210>
## 
## correctness rate = 0.8262
CVErr_LDA <- 1-unname(sat_LDA$result.pm["crossval.rate"])

Análisis Discriminante Cuadrático

set.seed(20)
sat_QDA <- stepclass(C ~ B1 + B2 + B3 + B4, data= sat_trn_df, method="qda", criterion = "CR", direction = "both", improvement <- 0.01, fold = 10)
##  `stepwise classification', using 10-fold cross-validated correctness rate of method qda'.
## 4435 observations of 4 variables in 6 classes; direction: both
## stop criterion: improvement less than 1%.
## correctness rate: 0.58107;  in: "B4";  variables (1): B4 
## correctness rate: 0.80496;  in: "B1";  variables (2): B4, B1 
## correctness rate: 0.84849;  in: "B2";  variables (3): B4, B1, B2 
## 
##  hr.elapsed min.elapsed sec.elapsed 
##        0.00        0.00        3.44
sat_QDA
## method      : qda 
## final model : C ~ B1 + B2 + B4
## <environment: 0x000001f8bde91e90>
## 
## correctness rate = 0.8485
CVErr_QDA <- 1-unname(sat_QDA$result.pm["crossval.rate"])

Comparación del Error de 10-fold CV

models <- c("tree", "rpart", "RF", "KNN", "RLM","LDA","QDA")
CVErrs <- c(CVErr_tree,CVErr_rpart,CVErr_rf,CVErr_KNN,CVErr_mlr,CVErr_LDA,CVErr_QDA)

CVErr_df <- data.frame(Label = models, Value = CVErrs)

ggplot(CVErr_df, aes(x = Label, y = Value)) + 
  geom_bar(stat = "identity", color = "black", fill = "red") + labs(x = "Modelos", y = "Error de Validación Cruzada K=10", title = "Comparación de Precisión de Modelos") +
  theme(plot.title = element_text(hjust = 0.5))

#conf_mat <- caret::confusionMatrix(pred_nnet, sat_tst_df$C)
#conf_mat
#conf_matrix(sat_tst_df$C,pred_nnet)
#probs <- predict(model_nnet, newdata = sat_tst_df, type = "probs")
#ROC <- multiclass.roc(response = sat_tst_df$C, predictor = probs)
#auc(ROC)

Evaluación de los Modelos contra el Set de Entrenamiento

Árbol de Clasificación

Librería Tree

Librería Rpart

Random Forest

Regresión Logística Multinomial

Vecinos Más Cercanos

Análisis Discriminante Lineal

Análisis Discriminante Cuadrático

ANEXO: Análisis Para el Set de Datos Completo

Queremos comprobar:

  • Que el píxel central es en efecto el principal predictor de la clase de la celda
  • Que despreciar los píxeles adyacentes al píxel central no produce una reducción significativa en el error de clasificación.
sat_trn_df_full <- data_trn
sat_trn_df_full$V37 <- as.factor(sat_trn_df_full$V37)

n <- nrow(sat_trn_df_full)
full_tree <- tree(V37 ~., data= sat_trn_df_full, control=tree.control(nobs= nrow(sat_trn_df_full), mincut = 1, minsize = 2))
full_tree_summary <- summary(full_tree)
EEnt <- full_tree_summary$misclass[1] / n
#full_tree_summary

set.seed(20)
full_tree_CV = cv.tree(full_tree, FUN = prune.misclass)
full_size <- full_tree_CV$size
full_dev <- full_tree_CV$dev
#ggplot(data=data.frame(full_size,full_dev), aes(x=full_size, y=full_dev)) +
#  geom_line(linetype="solid", color="blue", size=1.2)+
#  geom_point(color="red", size=3)+
#  labs(x= "Tamaño del Árbol", y = "Error de CV",title="Poda de Complejidad de Costos")

prune_full = prune.misclass(full_tree, best = 10)
plot(prune_full)
text(prune_full, pretty = 0)

set.seed(20)
full_tree_2 <-  rpart(V37 ~., data= sat_trn_df_full, method = "class", cp=-1,xval=10)
full_tree_2_summary <- printcp(full_tree_2)
## 
## Classification tree:
## rpart(formula = V37 ~ ., data = sat_trn_df_full, method = "class", 
##     cp = -1, xval = 10)
## 
## Variables actually used in tree construction:
##  [1] V1  V10 V11 V12 V13 V14 V15 V16 V17 V18 V19 V2  V20 V21 V22 V23 V24 V25 V26
## [20] V27 V28 V29 V3  V30 V31 V32 V33 V34 V35 V36 V4  V5  V6  V7  V8  V9 
## 
## Root node error: 3363/4435 = 0.75829
## 
## n= 4435 
## 
##             CP nsplit rel error  xerror      xstd
## 1   0.26167113      0   1.00000 1.00000 0.0084779
## 2   0.25334523      1   0.73833 0.76896 0.0097636
## 3   0.12905144      2   0.48498 0.48498 0.0095487
## 4   0.06363366      3   0.35593 0.35742 0.0088020
## 5   0.01754386      4   0.29230 0.29497 0.0082518
## 6   0.01486768      5   0.27475 0.28338 0.0081337
## 7   0.01159679      6   0.25989 0.25453 0.0078153
## 8   0.00683913      7   0.24829 0.24591 0.0077128
## 9   0.00669045      9   0.23461 0.23877 0.0076253
## 10  0.00624442     13   0.20458 0.23699 0.0076030
## 11  0.00490633     15   0.19209 0.21588 0.0073269
## 12  0.00356824     17   0.18228 0.21053 0.0072531
## 13  0.00297354     18   0.17871 0.20607 0.0071903
## 14  0.00267618     19   0.17574 0.20517 0.0071776
## 15  0.00237883     21   0.17038 0.20369 0.0071563
## 16  0.00223015     22   0.16800 0.19923 0.0070916
## 17  0.00208147     24   0.16354 0.19923 0.0070916
## 18  0.00178412     26   0.15938 0.19744 0.0070654
## 19  0.00168500     27   0.15760 0.19685 0.0070567
## 20  0.00163544     30   0.15254 0.19685 0.0070567
## 21  0.00148677     32   0.14927 0.19506 0.0070302
## 22  0.00133809     36   0.14332 0.19387 0.0070124
## 23  0.00118941     40   0.13797 0.19447 0.0070213
## 24  0.00104074     50   0.12608 0.19387 0.0070124
## 25  0.00089206     55   0.12073 0.19120 0.0069721
## 26  0.00077312     57   0.11894 0.19060 0.0069631
## 27  0.00074338     62   0.11508 0.19090 0.0069676
## 28  0.00059471     64   0.11359 0.19090 0.0069676
## 29  0.00049559     66   0.11240 0.19179 0.0069811
## 30  0.00029735     69   0.11091 0.19269 0.0069946
## 31  0.00019824     74   0.10943 0.19804 0.0070742
## 32  0.00014868     77   0.10883 0.19893 0.0070873
## 33  0.00000000     83   0.10794 0.19982 0.0071003
## 34 -1.00000000    157   0.10794 0.19982 0.0071003
full_tree_2_NP_pred <- predict(full_tree_2, sat_trn_df_full, type = "class")
# Poda
#full_tree_2$cptable 
#plotcp(full_tree_2) #grafico
#which.min(full_tree_2$cptable[,4])
cp<-full_tree_2$cptable[which.min(full_tree_2$cptable[,"xerror"]),"CP"]
poda<-prune(full_tree_2,cp=cp)
rpart.plot(poda, cex = 0.5, extra = 0)

set.seed(20)
rf_full<- randomForest(V37~., data= sat_trn_df_full, 
                     mtry=sqrt(ncol(sat_trn_df_full)-1), importance=TRUE)

#importance(rf_full)
varImpPlot(rf_full)

unname(rf_full$err.rate[nrow(rf_full$err.rate),1])
## [1] 0.08568207